{

VERSION 1.1, Daniel Backstr�m 2003, Ivan V. Minin 2008

http://www.baxtrom.com
http://code.google.com/p/matcalc

Using linlag.pas
----------------

The procedures defined in linalg.pas are designed to work on dynamic arrays
declared as "array of array of double". They can be accessed directly using
the square bracket notation (see the Delphi help for more info).

As of this version, the following functions and procedures are available:

   SetSize, NewMatrix, CopyMatrix, SizeOfMatrix, ZeroMatrix, Eye
   DelRow, DelColumn, SwapRows, SwapCols, ExtractSubMatrix, InsertSubMatrix
   Add, ScalarMultiply, Subtract,  Multiply
   Trace, MultTrace, Transpose, Diagonalize, Solve, Determinant
   SaveMatrixToFile, LoadMatrixFromFile

NOTE: Diagonalize, and thus indirectly Solve and Determinant requires that the
diagonal elements of the matrix are non-zero. It is up to the user to swap the
rows of the input matrice (and the right-hand-side vector in the Solve case) so
that its diagonal elements becomes non-zero.

more info on each function can be found below.


Some notes:
-----------

Solve returns nil if the system contains a zero diagonal element. To avoid this,
pivot the system. 


Example: solving a matrix system
--------------------------------

the system on augmented form (NOTE! Diagonal elements must be non-zero!):

| 3 -5 | 9 |
| 2  1 | 3 |

using linalg.pas in Delphi (presuming M and x are declared properly):

M := NewMatrix(SetSize(2, 2));
b := NewMatrix(SetSize(2, 1));

M[0, 0] := 3;
M[0, 1] := -5;
M[1, 0] := 2;
M[1, 1] := 1;

b[0, 0] := 9;
b[1, 0] := 3;

x := Solve(M, b);

..x should now be a 2x1 matrix containing x1 ~ 1.85 and x2 ~ -0.69


Copyright notice:
----------------

linalg, matrix algebra sourcecode for Delphi
Copyright (C) Daniel Backstr�m 2003

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

}
unit linalg;

interface
uses Dialogs;
type TMatrix = array of array of double;
type TSize = record
        Rows: Integer;
        Cols: Integer;
end;

function SetSize(rows, cols: Integer): TSize;
function NewMatrix(Size: TSize): TMatrix;
function CopyMatrix(M: TMatrix): TMatrix;
function SizeOfMatrix(M: TMatrix): TSize;
function ZeroMatrix(Size: TSize): TMatrix;
function Eye(Size: TSize): TMatrix;
function Add(M1, M2: TMatrix): TMatrix;
function ScalarMultiply(M: TMatrix; scalar: Double): TMatrix;
function Subtract(M1, M2: TMatrix): TMatrix;
function Multiply(M1, M2: TMatrix): TMatrix;

function DelRow(M: TMatrix; row: Integer): TMatrix;
function DelColumn(M: TMatrix; col: Integer): TMatrix;
function SwapRows(M: TMatrix; row1, row2: Integer): TMatrix;
function SwapCols(M: TMatrix; col1, col2: Integer): TMatrix;
function ExtractSubMatrix(M: TMatrix; StartRow, StartCol, NumRows, NumCols:
Integer): TMatrix;
function InsertSubMatrix(subM, M: TMatrix; StartRow, StartCol: Integer):
TMatrix;

function Trace(M: TMatrix): Double;
function MultTrace(M: TMatrix): Double;
function Transpose(M: TMatrix): TMatrix;
function Diagonalize(M: TMatrix): TMatrix;
function Solve(M, b: TMatrix): TMatrix;
function Determinant(M: TMatrix): Double;

procedure SaveMatrixToFile(M: TMatrix; FileName: String);
function LoadMatrixFromFile(FileName: String): TMatrix;

implementation


function int_min(a, b: Integer): Integer;
begin

if a < b then
        Result := a
else
        Result := b;

end;


function SetSize(rows, cols: Integer): TSize;
begin

Result.Rows := rows;
Result.Cols := cols;

end;


function NewMatrix(Size: TSize): TMatrix;
var
        i: Integer;
begin

SetLength(Result, Size.Rows);
for i := 0 to Size.Rows - 1 do
        SetLength(Result[i], Size.Cols);

end;


function CopyMatrix(M: TMatrix): TMatrix;
var
        Size: TSize;
        i, j: Integer;
begin

Size := SizeOfMatrix(M);
Result := NewMatrix(Size);

for i := 0 to Size.Rows - 1 do
        for j := 0 to Size.Cols - 1 do
        begin
                Result[i, j] := M[i, j];
        end;

end;


// Returns the size of a matrix
function SizeOfMatrix(M: TMatrix): TSize;
begin

SizeOfMatrix.Rows := length(M);
SizeOfMatrix.Cols := length(M[0]);

end;


// Returns zero matrix
function ZeroMatrix(Size: TSize): TMatrix;
var
        i, j: Integer;
begin

Result := NewMatrix(Size);

for i := 0 to Size.Rows - 1 do
        for j := 0 to Size.Cols - 1 do
                Result[i, j] := 0;

end;


// returns identity matrix
function Eye(Size: TSize): TMatrix;
var
        i, j: Integer;
begin

Result := NewMatrix(Size);

for i := 0 to Size.Rows - 1 do
        for j := 0 to Size.Cols - 1 do
        begin
                if i <> j then
                        Result[i, j] := 0
                else
                        Result[i, j] := 1;
        end;

end;


// Adding matrices
function Add(M1, M2: TMatrix): TMatrix;
var
        Size: TSize;
        i, j: Integer;
begin

Size := SizeOfMatrix(M1);

Result := NewMatrix(Size);

for i := 0 to Size.Rows - 1 do
        for j := 0 to Size.Cols - 1 do
                Result[i, j] := M1[i, j] + M2[i, j];

end;


// multiplication of scalar with matrix
function ScalarMultiply(M: TMatrix; scalar: Double): TMatrix;
var
        Size: TSize;
        i, j: Integer;
begin

Size := SizeOfMatrix(M);

Result := ZeroMatrix(Size);

if scalar = 0 then
        exit
else if scalar = 1 then
begin
        Result := CopyMatrix(M);
        exit;
end;

for i := 0 to Size.Rows - 1 do
        for j := 0 to Size.Cols - 1 do
                Result[i, j] := scalar*M[i, j];

end;


// matrix subtraction
function Subtract(M1, M2: TMatrix): TMatrix;
begin

Result := Add(M1, ScalarMultiply(M2, -1));

end;


// matrix multiplication
function Multiply(M1, M2: TMatrix): TMatrix;
var
        M1_Size, M2_Size, Res_Size: TSize;
        i, j, k: Integer;
begin

M1_Size := SizeOfMatrix(M1);
M2_Size := SizeOfMatrix(M2);
Res_Size.Rows := M1_Size.Rows;
Res_Size.Cols := M2_Size.Cols;

Result := ZeroMatrix(Res_Size);

for i := 0 to Res_Size.Rows - 1 do
        for j := 0 to Res_Size.Cols - 1 do
                for k := 0 to M1_Size.Cols - 1 do
                                Result[i, j] := Result[i, j] + M1[i, k]*
                                M2[k, j];

end;


// deletes the row:th row of M
function DelRow(M: TMatrix; row: Integer): TMatrix;
var
        Size, Res_Size: TSize;
        i, j: Integer;
begin

Size := SizeOfMatrix(M);
Res_Size.Rows := Size.Rows - 1;
Res_Size.Cols := Size.Cols;

Result := NewMatrix(Res_Size);

for i := 0 to Size.Rows - 1 do
        for j := 0 to Size.Cols - 1 do
        begin
                if i < row then
                        Result[i, j] := M[i, j]
                else if i > row then
                        Result[i - 1, j] := M[i, j];
        end;

end;


// deletes the col:th column of M
function DelColumn(M: TMatrix; col: Integer): TMatrix;
var
        Size, Res_Size: TSize;
        i, j: Integer;
begin

Size := SizeOfMatrix(M);
Res_Size.Rows := Size.Rows;
Res_Size.Cols := Size.Cols - 1;

Result := NewMatrix(Res_Size);

for i := 0 to Size.Rows - 1 do
        for j := 0 to Size.Cols - 1 do
        begin
                if j < col then
                        Result[i, j] := M[i, j]
                else if j > col then
                        Result[i, j - 1] := M[i, j];
        end;

end;


function SwapRows(M: TMatrix; row1, row2: Integer): TMatrix;
var
        Size: TSize;
        j: Integer;
        temp: Double;
begin

Result := CopyMatrix(M);

if row1 = row2 then Exit;

Size := SizeOfMatrix(M);

for j := 0 to Size.Cols - 1 do
begin
        temp := Result[row1, j];
        Result[row1, j] := Result[row2, j];
        Result[row2, j] := temp;
end;

end;


function SwapCols(M: TMatrix; col1, col2: Integer): TMatrix;
var
        Size: TSize;
        i: Integer;
        temp: Double;
begin

Result := CopyMatrix(M);

if col1 = col2 then Exit;

Size := SizeOfMatrix(M);

for i := 0 to Size.Rows - 1 do
begin
        temp := Result[i, col1];
        Result[i, col1] := Result[i, col2];
        Result[i, col2] := temp;
end;

end;


// Extracts submatrix from M
function ExtractSubMatrix(M: TMatrix; StartRow, StartCol, NumRows, NumCols:
Integer): TMatrix;
var
        Size: Tsize;
        i, j: Integer;
begin

Size := SizeOfMatrix(M);
Result := NewMatrix( SetSize(NumRows, NumCols) );

for i := StartRow to StartRow + NumRows - 1 do
        for j := StartCol to StartCol + NumCols - 1 do
                Result[i - StartRow, j - StartCol] := M[i, j];

end;


// Inserts submatrix into M
function InsertSubMatrix(subM, M: TMatrix; StartRow, StartCol: Integer):
TMatrix;
var
        Size: Tsize;
        i, j: Integer;
begin

Size := SizeOfMatrix(subM);

Result := M;
for i := StartRow to StartRow + Size.Rows - 1 do
        for j := StartCol to StartCol + Size.Cols - 1 do
                Result[i, j] := subM[i - StartRow, j - StartCol];

end;


// Returns trace (M11 + M22 + M33 + ..)
function Trace(M: TMatrix): Double;
var
        Size: TSize;
        i: Integer;
begin

size := SizeOfMatrix(M);

Result := 0;
for i := 0 to int_min(Size.Rows - 1, Size.Cols - 1) do
        Result := Result + M[i, i];

end;


// Returns multiplication trace (M11*M22*M33*..)
function MultTrace(M: TMatrix): Double;
var
        Size: TSize;
        i: Integer;
begin

size := SizeOfMatrix(M);

Result := 1;
for i := 0 to int_min(Size.Rows - 1, Size.Cols - 1) do
        Result := Result*M[i, i];

end;


// Returns the transpose of M
function Transpose(M: TMatrix): TMatrix;
var
        Size: TSize;
        i, j: Integer;
begin

Size := SizeOfMatrix(M);

// Swapping size
i := Size.Rows;
Size.Rows := Size.Cols;
Size.Cols := i;

Result := NewMatrix(Size);

for i := 0 to Size.Rows - 1 do
        for j := 0 to Size.Cols - 1 do
                Result[i, j] := M[j, i];

end;


{ Returns the diagonalized form of M (can be non-square, i.e. an augmented
system }
function Diagonalize(M: TMatrix): TMatrix;
var
        Size: TSize;
        i, j: Integer;
        Row, RefRow: TMatrix;
        alpha: Double;
begin

// returns nil if any of the diagonal elements is zero
if MultTrace(M) = 0 then
begin
        Result := nil;
        Exit;
end;

Size := SizeOfMatrix(M);

Result := CopyMatrix(M);

for i := 0 to Size.Rows - 1 do
begin
        RefRow := ExtractSubMatrix(Result, i, 0, 1, Size.Cols);

        for j := 0 to Size.Rows - 1 do
        begin
                if i <> j then
                begin
                        Row := ExtractSubMatrix(Result, j, 0, 1, Size.Cols);

                        alpha := -Row[0, i]/RefRow[0, i];
                        Row := Add( Row, ScalarMultiply(RefRow, alpha) );

                        Result := InsertSubMatrix(Row, Result, j, 0);
                end;
        end;
end;

// if augmented system (not square)
if Size.Rows <> Size.Cols then
        for i := 0 to Size.Rows - 1 do
        begin
                Result[i, Size.Cols - 1] := Result[i, Size.Cols - 1]/
                Result[i, i];
                Result[i, i] := 1;
        end;

end;


// Solves system Mx = b using function Diagonalize (diagonal elements of M <>0)
function Solve(M, b: TMatrix): TMatrix;
var
        Size: TSize;
        AugmentedSystem: TMatrix;
begin

Size := SizeOfMatrix(M);
Size.Cols := Size.Cols + 1;
AugmentedSystem := NewMatrix( Size );
Result := NewMatrix( SizeOfMatrix(b) );

InsertSubMatrix(M, AugmentedSystem, 0, 0);
InsertSubMatrix(b, AugmentedSystem, 0, Size.Cols - 1);

AugmentedSystem := Diagonalize(AugmentedSystem);
Result := ExtractSubMatrix( AugmentedSystem, 0, Size.Cols - 1, Size.Rows, 1);

end;


// Returns the determinant of M
function Determinant(M: TMatrix): Double;
begin

Result := MultTrace( Diagonalize( M ) )

end;


procedure SaveMatrixToFile(M: TMatrix; FileName: String);
var
        Size: TSize;
        Rows, Cols: Double;
        i, j: Integer;
        F: file of double;
        sign: Double;
begin

Size := SizeOfMatrix(M);

AssignFile(F, FileName);
Rewrite(F);

Rows := Size.Rows;
Cols := Size.Cols;
sign:=14613180;

Write(F, sign);
Write(F, Rows);
Write(F, Cols);

for i := 0 to Size.Rows - 1 do
        for j := 0 to Size.Cols - 1 do
                Write(F, M[i, j]);

CloseFile(F);

end;


function LoadMatrixFromFile(FileName: String): TMatrix;
var
        Size: TSize;
        Rows, Cols, temp: Double;
        i, j: Integer;
        F: file of double;
        sign: Double;
begin

AssignFile(F, FileName);
Reset(F);

Read(F, sign);
if sign<>14613180 then begin Dialogs.ShowMessage('���� ���������!!!'); exit; end else begin
Read(F, Rows);
Read(F, Cols);

Size.Rows := Round( Rows );
Size.Cols := Round( Cols );

Result := NewMatrix( Size );

for i := 0 to Size.Rows - 1 do
        for j := 0 to Size.Cols - 1 do
        begin
                Read(F, temp);
                Result[i, j] := temp;
        end;

CloseFile(F);
end;
end;

end.
